\documentclass[11pt,a4paper]{book}
\usepackage{amsmath,amssymb}
\usepackage[utf8]{inputenc}
\usepackage[french]{babel}
\usepackage[pdftex]{graphicx}
\usepackage[T1]{fontenc}
\usepackage[top=2cm, bottom=2cm, left=2cm, right=2cm]{geometry}
\usepackage{fancyhdr}
\usepackage{wasysym}

\def\blurb{%
 Institut National Polytechnique de Toulouse \\ 
 ENSEEIHT \\
 Département d'Informatique et Mathématiques Appliquées \\[1em]
 Rapport de projet d'équations aux dérivées partielles}
\def\clap#1{\hbox to 0pt{\hss #1\hss}}%
\def\ligne#1{%
  \hbox to \hsize{%
    \vbox{\centering #1}}}%
\def\haut#1#2#3{%
  \hbox to \hsize{%
    \rlap{\vtop{\raggedright #1}}%
    \hss
    \clap{\vtop{\centering #2}}%
    \hss
    \llap{\vtop{\raggedleft #3}}}}%
\def\bas#1#2#3{%
  \hbox to \hsize{%
    \rlap{\vbox{\raggedright #1}}%
    \hss
    \clap{\vbox{\centering #2}}%
    \hss
    \llap{\vbox{\raggedleft #3}}}}%

\begin{document}

\thispagestyle{empty}\vbox to 1\vsize{%
  \vss
  \vbox to 1\vsize{%
    \haut{}{\blurb}{}
    \vfill
    \vspace{1cm}
    \begin{flushleft}
	\Huge Etude d'un problème de diffusion évolutif
    \end{flushleft}
    \par
    \hrule height 2pt
    \par
    \begin{flushright}
	\Large Yuravin Bung, Hugo Charton
	\par
    \end{flushright}
    \vspace{1cm}
    \vfill \vfill
    \ligne{%
      \hfill
      \begin{tabular}{l}
        Avril 2011
      \end{tabular}
      }
    }%
  \vss
  }

\pagestyle{fancy}
\renewcommand{\chaptermark}[1]{\markboth{#1}{}}
\renewcommand{\sectionmark}[1]{\markright{\thesection\ #1}}
\fancyhf{} \fancyhead[LE,RO]{\bfseries\thepage}
\fancyhead[LO]{\bfseries\rightmark}
\fancyhead[RE]{\bfseries\leftmark}
\renewcommand{\headrulewidth}{0.5pt}
\addtolength{\headheight}{0.5pt}
\renewcommand{\footrulewidth}{0pt}
\fancypagestyle{plain}{ \fancyhead{}
\renewcommand{\headrulewidth}{0pt}}

\newpage $\text{ }$
\renewcommand{\headrulewidth}{0.4pt}
\fancyhead[LO]{\bfseries\rightmark}
\fancyhead[RE]{\bfseries\leftmark}
\tableofcontents 


\chapter{Introduction}

\section{Position du problème}

Le but de ce projet est de résoudre le problème de diffusion évolutif suivant :

$$
\left\{
\begin{array}{r c l}
\cfrac{{\partial}u(x,y,t)}{{\partial}t} - {\nu}{\Delta}u(x,y,t) + a\cfrac{{\partial}u(x,y,t)}{{\partial}x} + b\cfrac{{\partial}u(x,y,t)}{{\partial}y} = f(x,y,t)\\
u(x,y,t)_{|{\partial}\Omega} = 0 : \emph{conditions limites} \\
u(x,y,0) = u_0(x,y) : \emph{conditions initiales}
\end{array}
\right.
$$

\section{Axes d'études}

Pour faciliter la résolution de ce problème, il était conseillé de le découper en deux blocs :
\begin{itemize}
\item L'étude dans un premier temps du problème de diffusion stationnaire : sans la partie temporelle.
\item Dans un deuxième temps du prolème de diffusion évolutif.
\end{itemize}

En effet, si ce découpage n'avait pas été fait, il aurait été très chanceux de notre part de réussir à débugger, puis à faire marcher notre projet en temps voulu !


\chapter{Méthode de Jacobi, méthode de Gauss-Seidel}

\section{Introduction}

Nous allons donc d'abord étudier la partie stationnaire avec les deux méthodes :
\begin{itemize}
\item Jacobi
\item Gauss-Seidel
\end{itemize}
Le problème se ramène à celui-ci :
$$
\left\{
\begin{array}{r c l}
- {\epsilon}{\Delta}u(x,y) + {\alpha}\cfrac{{\partial}u(x,y,t)}{{\partial}x} + {\beta}\cfrac{{\partial}u(x,y,t)}{{\partial}y} + cu(x,y) = g(x,y)\\
u(x,y,t)_{|{\partial}\Omega} = 0 : \emph{conditions limites} \\
\end{array}
\right.
$$

Cela revient donc à résoudre l'équation :

$Au = g$

Suite à une discrétisation faite, nous avons constaté que la matrice avait ce format :
\[
A = 
\begin{pmatrix}

   \begin{pmatrix}
	\cfrac{4{\epsilon}}{h^2} + c & \cfrac{\beta}{2h} - \cfrac{\epsilon}{h^2} & \cdots \\
	-\cfrac{\epsilon}{h^2} - \cfrac{\beta}{2h} & \ddots &\vdots \\
	\vdots & \cdots & \vdots 
	\end{pmatrix}

 & \begin{pmatrix}
    -\cfrac{\epsilon}{h^2} + \cfrac{\alpha}{2h} & \cdots \\
    \vdots & \cdots & \vdots 
	\end{pmatrix}
 & \cdots \\
 
	\begin{pmatrix}
    -\cfrac{\epsilon}{h^2} - \cfrac{\alpha}{2h} & \cdots \\
    \vdots & \cdots & \vdots 
	\end{pmatrix}
	& \ddots \\
	
   \vdots & \ddots \\
\end{pmatrix}
\]

En choisissant $u$ de la forme :
\[
u = 
\begin{pmatrix}
	u_{11} \\
	u_{12} \\
	\vdots \\
	u_{1n} \\
	u_{21} \\
	\vdots 
\end{pmatrix}
\]

A est bien sûr constituée de $n$ blocs de taille $n$.
Par soucis d'espace mémoire, nous avons donc préféré stocker la matrice A sous cette forme, qui est beaucoup plus légère au niveau de la mémoire :
\[
A = 
\begin{pmatrix}
	-\cfrac{\epsilon}{h^2} + \cfrac{\alpha}{2h} \\
	\cfrac{\beta}{2h} - \cfrac{\epsilon}{h^2} \\
	\cfrac{4{\epsilon}}{h^2} + c \\
	-\cfrac{\epsilon}{h^2} - \cfrac{\beta}{2h} \\
	-\cfrac{\epsilon}{h^2} - \cfrac{\alpha}{2h}
\end{pmatrix}
\]

Il a donc fallu adapter le \textbf{produit matrice/vecteur} en fonction de la forme de la matrice $A$.

La solution exacte de référence est : 
\newline
$ u_{ex}(x,y) = xy(x-1)(y-1)$
\newline
Et en remplaçant dans le système, on obtient ainsi le second membre suivant :
\newline
$g(x,y) = -2{\epsilon}(y(y-1)+x(x-1)) + {\alpha}y(y-1)(2x-1) + {\beta}x(x-1)(2y-1) + cxy(x-1)(y-1)$

Trois différents jeux de données pour tester les deux méthodes précédemment citées sont demandées d'utiliser :
\begin{enumerate}
\item L'équation de diffusion.
\item L'équation de réaction-convection-diffusion.
\item L'équation de convection-diffusion avec convection prépondérante.
\end{enumerate}

Nous allons maintenant voir les résultats de chaque méthode.

\section{Schéma de Jacobi}

La méthode de \textbf{Jacobi} se base donc sur cette équation :
\newline
\textbf{$u_{n+1} = D^{-1}(g - Ru_n)$}
avec 
\begin{itemize}
\item $D$ la diagonale de $A$ (donc, ici, un scalaire) 
\item et $R$ le reste de la matrice $A$ (un vecteur de dimension 4).
\end{itemize}


Voici les différents résultats obtenus :
\newline
\begin{tabular}{|c|c|c|c|}
\hline
$n=100$, $seuilResidu = 10^{-6}$ & Nombre d'itérations & Temps de résolution & Erreur \\
\hline
Diffusion & 25469 & 12.19s & 9.9989E-007\\
\hline
Réaction-convection-diffusion & 20373 & 9.87s & 9.9949E-007 \\
\hline
Convection(prépondérante)-diffusion &  339 &  0s & Nan \\
\hline
\end{tabular}

\section{Schéma de Gauss-Seidel}

Alors que la méthode de \textbf{Gauss-Seidel} se base sur cette équation :
\newline
\textbf{$u_{n+1} = L^{-1}(g - Hu_n)$}
avec 
\begin{itemize}
\item $L$ la matrice triangulaire inférieure contenant la diagonale (donc vecteur de dimension 3) 
\item et $H$ la matrice triangulaire supérieure (donc un vecteur de taille 2)
\end{itemize}


Voici les différents résultats obtenus :
\newline
\begin{tabular}{|c|c|c|c|}
\hline
$n=100$, $seuilResidu = 10^{-6}$ & Nombre d'itérations & Temps de résolution & Erreur \\
\hline
Diffusion & 26870 & 10.82s & 9.9975E-007\\
\hline
Réaction-convection-diffusion & 18593 & 7.51s & 9.9946E-007 \\
\hline
Convection(prépondérante)-diffusion &  9 &  0s & Nan \\
\hline
\end{tabular}

On constate donc que la méthode Gauss-Seidel est légèrement plus rapide que la méthode de Jacobi, alors qu'elles font à peu près le même nombre d'itérations.
Nous verrons dans la suite (résolution du problème évolutif) qu'il est alors plus intéressant d'utiliser la méthode de Gauss-Seidel.


\chapter{Schémas de Euler, Crank-Nicholson}

\section{Introduction}

Après avoir étudié la résolution du problème de diffusion stationnaire, nous allons maintenant étudier le problème général.
Plusieurs schémas sont étudiés :
\begin{enumerate}
\item \textbf{Méthode d'Euler explicite} qui consite à calculer :
\newline
$ \cfrac{u_{n+1} - u_n}{{\Delta}t} = h(t_n,u_n)$
\item \textbf{Méthode de Crank-Nicholson} qui consiste à calculer :
\newline
$ \cfrac{u_{n+1} - u_n}{{\Delta}t} = \cfrac{1}{2}(h(t_n,u_n)+h(t_{n+1},u_{n+1}))$
\end{enumerate}

La solution exacte du problème est donné : 
\newline
$u_{ex}(x,y,t) = \exp{(-{\pi}t)}\sin{({\pi}x)}\sin{({\pi}y)}$

Grâce à cette solution exacte, il est donc possible de calculer le second membre $f$ !
En effet, il suffit de remplacer cette fonction dans l'équation \textbf{(1)} du sujet, de calculer les différentes dérivées, selon le temps, le déplacement x, puis y et enfin de Laplacien.
\newline
On obtient ainsi :
$f(x,y,t) = {\pi}{\exp{(-{\pi}t)}}(-{\sin{({\pi}x})}{\sin{({\pi}y})} + a{\cos{({\pi}x})}{\sin{({\pi}y})} + b{\sin{({\pi}x})}{\cos{({\pi}y})} - 2{\nu}{\pi}{\sin{({\pi}x})}{\sin{({\pi}y})})$


Maintenant que la fonction second membre $f$ est obtenue, le code précèdemment expliqué peut être utilisé pour inspirer ceux qui vont arriver.
Seul un détail a été modifié, en effet, nous avions implémenté une \emph{subroutine} qui créait le second membre $g$, cependant pour éviter de créer une matrice second membre, nous avons inséré ce code dans une boucle temporelle.
Tout va être expliqué plus en détail dans ce qui suit.

\section{Schéma d'Euler}

Comme on vient de le dire, \textbf{la méthode d'Euler explicite} calcule : 
\newline
$ \cfrac{u_{n+1} - u_n}{{\Delta}t} = h(t_n,u_n)$
\newline
On obtient donc :
$u_{n+1} = {\Delta}t(f-Au_n) + u_n$
\newline
Avec A la matrice habituelle vue dans la premier partie, seules ses coefficients changent :
\[
A = 
\begin{pmatrix}
\cfrac{a}{2h} - \cfrac{\nu}{h^2} \\
\cfrac{b}{2h} - \cfrac{\nu}{h^2} \\
\cfrac{4{\nu}}{h^2} \\
- \cfrac{b}{2h} - \cfrac{\nu}{h^2} \\
- \cfrac{a}{2h} - \cfrac{\nu}{h^2}
\end{pmatrix}
\]
Cette méthode est temporelle, et nous voulons étudier plusieurs valeurs selon un ensemble de temps assez élevé.
En effet, si le découpage temporel (temps de discrétisation) n'est pas assez élevé, la convergence n'est pas atteinte.

Par exemple, avec deux itérations temporelles, pour une matrice de trois blocs de taille trois, on obtient pour \textbf{l'équation de diffusion} :
\begin{itemize}
\item Nombre d'itérations : 10000 (la limite que nous avons imposée).
\item Valeur du résidu : 2.9864789866893995
\end{itemize}

Il est normal que le nombre d'itérations atteigne la limite puisque la convergence n'est pas atteinte.

Cependant, nous n'avons pas réussi à coder cette partie sans écrire \textbf{write(*,*) u}, et le temps de résolution se trouve vite handicapé, puisqu'il doit afficher $u$ dans chaque pas d'itération.

\section{Schéma de Crank-Nicholson}

Contrairement à la méthode d'Euler, \textbf{la méthode de Crank-Nicholson implicite}  utilise cette discrétisation :
\newline
$u_{n+1} = {\Delta}t(f-Au_n) + u_n$
\newline
Et on obtient ainsi :
\newline
$(I + \cfrac{{\Delta}t}{2}A)u_{n+1} = \cfrac{{\Delta}t}{2}(f_n - Au_n + f_{n+1}) + u_n$
\newline
Avec A la même matrice que pour la méthode d'Euler, et I la matrice identité de $n$ blocs de taille $n$.
Cependant, vu la forme de la matrice A, cela revient à ajouter $1$ au terme central de $A$ : $\cfrac{4{\nu}}{h^2}$.

On contaste cependant que cette fois, l'équation obtenue est \emph{explicite}, c'est-à-dire que $u_{n+1}$, que l'on doit déterminer ne peut pas être trouvé directement.
En effet, il est d'abord multiplié par une matrice, c'est pourquoi, afin de résoudre ce problème, il faut utiliser les méthodes stationnaires vues dans la première partie :
\begin{itemize}
\item la méthode de Jacobi
\item celle de Gauss-Seidel
\end{itemize}

De plus, comme la deuxième méthode citée est plus efficace (en temps de résolution, pas en nombre d'itérations) que la première, nous avons trouvé judicieux d'utiliser plutôt la méthode de \textbf{Gauss-Seidel} dans cette résolution.


Nous avons constaté le même problème que pour la méthode d'Euler, au niveau de la résolution, et n'avons donc pas pu obtenir de résultats, même si au niveau algorithme, nous avons compris.


\chapter{Conclusion}

Ce projet était très intéressant, en effet, nous avons du résoudre, d'abord les équations de la chaleur stationnaire, puis le problème évolutif.
Cependant, faute de temps, nous n'avons pas pu aboutir la partie évolutive.

En effet, nous avons perdu beaucoup de temps dans le débuggage du code en Fortran.
C'est un des principal problème de ce projet : le langage de programmation.
Celui-ci est assez capricieux, par exemple, le résultat donné peut être Infinity, alors que si un \textbf{write(*,*) ''} judicieusement ajouté permet de donner un résultat cohérent !
Il est assez frustrant de ne pas pouvoir terminer le projet à cause de petits détails comme celui-ci.


\end{document}
